/***********************************************************************
               _                     
              | |                    
 _ __    __ _ | | _ __ ___       ___ 
| '_ \  / _` || || '_ ` _ \     / __|
| |_) || (_| || || | | | | | _ | (__ 
| .__/  \__,_||_||_| |_| |_|(_) \___|
| |                                  
|_|                                  

palm.c
Mark Bolding
and a cast of thousands
Yun Jiang
D Twieg
R Deshpande
Maj Hedehus
etc.

adapted from Rose.c adapted from episs.c
waveforms are now pre-computed externally in matlab or your favorite enviroment
 *
 * revised 6 dec 2006 -- msb 
 * remove more unused variables
 * clean up the formatting
 * add diagnostic ouput

 * revised from pinecone.c 29 oct 2007
 * add prespin echo gradient and multishot
 *


*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
     _                       __                                 __     _
    |                       |                                     |     |
RF  | -----[90deg]----------|---------------- |** ACQ **|---------|     |
    |                       |                                     |     |
    |        ____           |                                     |     |
GSS | ______/////\     _____|_____________________________________|     |
    | <---->      \\\\/ gssr|                                     |     |
    |   predelay  <--> tref |                                     | etl | ns
    |                       |                                     |     |
    |                       |                                     |     |
GRO |  _____________________|_________________//////////\_________|     |
    |                       |                                     |     |
    |                       |                                     |     |
    |                       |                                     |     |
GRO |  _____________________|_________________//////////\_________|     |
    |                       |                                     |     |
    |_                      |__                                 __|    _|

*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/

#include "standard.h"
#include "group.h"
#define MAX_DAC_VAL  32767.0  /* maximum DAC value for system*/

pulsesequence()
{
    /* INTERNAL VARIABLE DECLARATIONS *********************/
	double	sefudge;
	double  sedelay,predelay,seqtime;
	double  agss,grate,gssrint;
	double  rgss,rgssr;
	double  ted,tssr;
	double  agssr,tweakspoil;
	double  grad_dur,grad_pre_dur; /* length of the shaped gradients in seconds*/
	char    ro_shape[MAXSTR], pe_shape[MAXSTR];       /* base name of gradient tables */
	char    grad_ro[MAXSTR], grad_pe[MAXSTR];         /* full path of grad tables */
	char    ro_pre_shape[MAXSTR], pe_pre_shape[MAXSTR]; /* fid (preSE 180deg excite) grad tables */
	char    grad_pre_ro[MAXSTR], grad_pre_pe[MAXSTR];
	char    *home_dir, calib[MAXSTR], spinecho[MAXSTR], fatsat[MAXSTR];

	sefudge=getval("sefudge");
    
	initparms_sis();
	getstr("fatsat", fatsat);     /*fat saturation*/ 
	getstr("calib",   calib);     /* for slice selctive cal with big phantoms */ 
	getstr("spinecho",spinecho);  /* spin echo/gradient echo flag */
	grad_dur = getval("grad_dur");/* get the shaped gradient durations */
	grad_pre_dur = getval("grad_pre_dur");
    tweakspoil = getval("tweakspoil");
	grate = trise/gmax;           /* gradient rate */
	agss = fabs(gss);             /* absolute val of slice sel gradient */
	rgss = grate*agss;            /* slice select gradient rise time */

  	/* GRADIENT*TIME INTEGRALS FOR SLICE SELECTION, PHASE AND READOUT */
	gssrint = agss*(p1/2.0 + rof2 + rgss/2.0); /* refocus integral */
	gssr = gss*gssf;      /* gssf: slice select refocus fraction */
	agssr = fabs(gssr);
	rgssr = grate*agssr;
	tssr = gssrint/agssr; /* slice refocussing time for constant pulse */
	tssr = tssr - rgssr;  /* correct for ramps */
    
    /* load the gradient tables */
	if ((home_dir = (char *) getenv("HOME")) == NULL) 
	{
		abort_message("Error in determining HOME environment variable\n");
	}
	getstr("ro_shape",ro_shape);
	getstr("pe_shape",pe_shape);
	getstr("ro_pre_shape",ro_pre_shape);
	getstr("pe_pre_shape",pe_pre_shape);

	sprintf(grad_ro, "%s.GRD",ro_shape);
	sprintf(grad_pe, "%s.GRD",pe_shape);
	sprintf(grad_pre_ro, "%s.GRD",ro_pre_shape);
	sprintf(grad_pre_pe, "%s.GRD",pe_pre_shape);

  	/* time bookkeeping */
	seqtime = rgss*2.0 + rof1 + p1 + rof2 + rgssr*2.0 + tssr + np/(2*sw); /* rof: reciever turn on delays */
    ted = 0.0; /* echo delay */
	if (spinecho[0] == 'y') /* add 180 SE excitation to seq duration */
	{	
		sedelay = rof1 + p2 + rof2 + (rgss*2.0);
        ted = p1/2.0 + rgss + rof2 + 2.0*rgssr + tssr + np/(6.0*sw) + rgss + tspoil + p2/2.0;
		seqtime = seqtime + sedelay;
	}
    if (te/2 - ted < 0.0)
    {
        abort_message("%s: Requested te too short. Min te = %f ms",seqfil,2.*ted*1000.);
    }
	
    if (ir[0] == 'y') /* add 180 IR excitation to seq duration */
	{		
		sedelay = rof1 + p2 + rof2 + (rgss*2) + (tspoil*2);
		seqtime = seqtime + sedelay;
	}
    
	if (fatsat[0] == 'y') /* add in fat saturation time */
		seqtime = seqtime + rof1 + psat + rof2 + tcrush + trise; 		
	
	if (tr - seqtime < 0.0) 
	{
		abort_message("%s: Requested tr too short. Min tr = %f ms",seqfil,seqtime*1000.);
	}
	predelay = (tr - seqtime);
    
	/* if (verbose[0] == 'y') 
	{ /* print test information */
		printf("*********** %s ************************************\n",seqfil);
		printf("gradient files\n ro: %s\n pe: %s\n", grad_ro, grad_pe);
		printf("gradient files\n se ro: %s\n se pe: %s\n", grad_pre_ro, grad_pre_pe);
		printf("calibration setting: %s\n", calib);
		printf("gro: %.3e\n", gro);
		printf("gpe: %.3e\n", gpe);
		printf("gss: %.3e\n", gss);
		printf("gssr: %.3e\n", gssr);
		printf("at: %.9f\n", at);
		printf("ted: %.9f\n", ted);
        printf("te: %.9f\n" , te);
        printf("p1: %.9f\n" , p1);
        printf("p2: %.9f\n" , p2);
        printf("rof1: %.9f\n" , rof1);
        printf("rof2: %.9f\n" , rof2);
        printf("np/(2*sw): %.9f\n" , np/(2*sw));
        printf("np: %.9f\n" , np);
        printf("np/3: %.9f\n" , np/3.);
        printf("2*np/3: %.9f\n" , 2*np/3);        
		printf("***************************************************\n");
    /*}*/
  
  	/* PULSE SEQUENCE *************************************/
    /* Relaxation delay ***********************************/       
    status(A);
    /* Phase cycle: Part I ********************************/
    mod2(ct,v3);                    /*0101 0101*/
    dbl(v3,v3);                     /*0202 0202*/
    hlv(ct,v4);                     /*0011 2233*/
    mod2(v4,v4);                    /*0011 0011*/
    add(v3,v4,v1);                  /*0213 0213*/
    assign(v1,oph);
    add(one,v1,v2);
    assign(zero,v10);
    assign(zero,v12);
    
	delay(predelay);       
	msloop(seqcon[1],ns,v13,v14); /* Slice loop ***************/       
		rotate();              
		/* fat saturation pulse *******************************/
		if (fatsat[0] == 'y') 
		{
			obspower(satpwr);
			obsoffset(resto-(3.4*sfrq)); /* hardcoded fat sat offset */
			shapedpulse(satpat,psat,zero,rof1,rof2);
			obl_gradient(gcrush,gcrush,gcrush);
			delay(tcrush);
			zero_all_gradients();
			delay(trise);
			obsoffset(resto);
		}      
        
		/* Gating *********************/
		/* xgate(ticks); */
        
		poffset_list(pss,gss,ns,v14);  /* set slice offset frequency */	
		if (ir[0] == 'y')
		{
			/* inversion 180 degree RF pulse ******************/
			obspower(tpwr2);
			obl_gradient(0.0,0.0,gss);
			delay(rgss); /* rise time */
			delay(tspoil);
			shapedpulse(p2pat,p2,v2,rof1,rof2);
			delay(tspoil);
			zero_all_gradients();
			delay(rgss);		
			delay(ti-tspoil-rgss-p2/2); /* wait for the inversion time */
		}
		/* 90 degree RF excitation pulse ****************************/ 
		observepower(tpwr1);
		if (calib[0] == 'r') obl_gradient(gss,0.0,0.0);
		if (calib[0] == 'p') obl_gradient(0.0,gss,0.0);
		if (calib[0] == 'n') obl_gradient(0.0,0.0,gss);
		delay(rgss);
		shapedpulse(p1pat,p1,v1,rof1,rof2); /* actual excitation */
		zero_all_gradients();
		delay(rgss);	
		
		/* Slice refocus  ********/
		if (calib[0] == 'r') obl_gradient(-gssr,0.0,0.0);
		if (calib[0] == 'p') obl_gradient(0.0,-gssr,0.0);
		if (calib[0] == 'n') obl_gradient(0.0,0.0,-gssr);
		delay(rgssr);
		delay(tssr); /* refocus time */
		zero_all_gradients();
		delay(rgssr);
		
		/* spin echo 180 degree RF pulse *****************************/
		if (spinecho[0] == 'y')
		{

			if (calib[0] == 'r') obl_shapedgradient(grad_pre_ro,"","",grad_pre_dur,gro,0.0,0.0,1,NOWAIT);
			if (calib[0] == 'p') obl_shapedgradient("",grad_pre_pe,"",grad_pre_dur,0.0,gpe,0.0,1,NOWAIT);
			if (calib[0] == 'n') 
			{
				obl_shapedgradient(grad_pre_ro,grad_pre_pe,"",grad_pre_dur,gro,gpe,0.0,1,NOWAIT);
	        }
			/* delay(ted); /* delay is taken up by acquire? */
			obsoffset(resto); /* put the reciever back */
			acquire(np/3.,1./sw); /* acquire data from fid before spin echo */
			delay(te/2-ted);
			obspower(tpwr2);
            poffset_list(pss,gss,ns,v14); /* move transmitter to slice offset frequency */
            if (calib[0] == 'r') obl_gradient(1*gss,0.0,0.0);
			if (calib[0] == 'p') obl_gradient(0.0,1*gss,0.0);
			if (calib[0] == 'n') obl_gradient(0.0,0.0,1*gss);
            delay(rgss);
            delay(tspoil);
			if (calib[0] == 'r') obl_gradient(gss,0.0,0.0);
			if (calib[0] == 'p') obl_gradient(0.0,gss,0.0);
			if (calib[0] == 'n') obl_gradient(0.0,0.0,gss);
			shapedpulse(p2pat,p2,v2,rof1,rof2); /* 180 excitation */
            if (calib[0] == 'r') obl_gradient(1*gss,0.0,0.0);
			if (calib[0] == 'p') obl_gradient(0.0,1*gss,0.0);
			if (calib[0] == 'n') obl_gradient(0.0,0.0,1*gss);
            delay(tspoil*tweakspoil);
			zero_all_gradients();
			delay(rgss);
		}
		
		if (calib[0] == 'r') obl_shapedgradient(grad_ro,"","",grad_dur,gro,0.0,0.0,1,NOWAIT);
		if (calib[0] == 'p') obl_shapedgradient("",grad_pe,"",grad_dur,0.0,gpe,0.0,1,NOWAIT);
		if (calib[0] == 'n') 
        {
			obl_shapedgradient(grad_ro,grad_pe,"",grad_dur,gro,gpe,0.0,1,NOWAIT);
        }
		obsoffset(resto);         /* put the reciever back */
        if (spinecho[0] == 'y')   /* show me the data! */
        {
            acquire((sefudge+2.*np/3.),1./sw);    
        }
        else
            acquire(np,1./sw);
		zero_all_gradients();
	endmsloop(seqcon[1],v14);
}
